/*
    本实验根据数组x[], y[]列出的一组数据，用最小二乘法求它的拟合曲线。
    近似解析表达式为y = a0 + a1 * x + a2 * x^2 + a3 * x^3;
*/
#include <stdio.h>
#include <math.h>

#define maxn 100
#define rank_ 3


    static double calc(double x, double b0, double b1, double b2, double b3)
    {
        return (b0+ b1*x +b2*x*x+b3*x*x*x);
    }

    static double Calc(double start, double end, double b0, double b1, double b2, double b3)
    {
        double sum = 0;
        const int N = 100000;
        const double dx = (end - start) / N;
        for (int i = 0; i < N; ++i)
        {
            sum += (calc(start + dx * i, b0,b1,b2,b3) * dx);
        }
        return sum;
    }

int main(){
    double x[] = { 622.0779755,610.370549,600.2421832,585.4214305,567.726079,557.5426511,544.1358074,528.7464428,518.2136065,501.4755537,486.8921707,473.3578461,462.2417877,445.0489251,433.3557722,420.309919,407.4950751,391.1955973,380.7625003,363.7512639,352.8362616,337.8004148,324.260974,309.3247626,295.9831238,281.1527377,269.034056,254.1283066,240.7785149,226.7522598,214.4607141,200.9722862,187.1517744,172.7814125,160.1647132,145.1062564,131.692064,117.3713136,103.8473297,91.50472996,77.51463196,63.47121748,50.26840197,36.4121351,22.06140997,8.721003824,-5.333364023,-19.40743067,-33.10927729,-47.19096787,-60.98378011,-74.57759345,-88.92301317,-102.6209973,-115.6427795,-131.1713465,-145.0562804,-159.1138717,-173.1656689,-186.4281611,-200.4131584,-212.9764586,-224.89985,-238.4668604,-254.1886764,-269.0459643,-282.6854897,-294.362019,-308.8986245,-321.4411213,-336.5133568,-350.4103157,-362.7218684,-377.9563789,-392.6527258,-405.9347226,-418.8557091,-434.1725962,-449.5604938,-463.7019719,-476.0180185,-492.1879172,-501.3179633,-520.6196803,-535.4507721,-549.5302043,-557.6047441,-573.1328557,-585.4059672,-600.935568,-613.961953,-626.9254239,-642.526168,-654.5085784,-670.126434,-684.3542987,-698.6106628,-714.5664736,-702.1294032,-716.918282 };
    double y[] = { 1260.894521,1268.822995,1281.379851,1281.625042,1277.533209,1289.639947,1293.184141,1293.072387,1304.86423,1299.604274,1300.867793,1304.789772,1314.070595,1307.322628,1316.510074,1321.782347,1327.881306,1322.349048,1335.792244,1326.012827,1338.268498,1335.95467,1340.33422,1337.702206,1341.741402,1338.796899,1346.383926,1346.223906,1351.722866,1352.118491,1360.196898,1364.276783,1366.241272,1360.068963,1371.680817,1362.293718,1367.674377,1357.937029,1362.046891,1384.980464,1377.821281,1376.537469,1384.087457,1390.523339,1373.822876,1387.972602,1388.989761,1389.864508,1394.607033,1393.200995,1385.65868,1385.995015,1390.158875,1384.201189,1377.153131,1392.837061,1389.448695,1387.909138,1390.257045,1382.486724,1380.528727,1372.574963,1358.509867,1355.178791,1363.509119,1368.809435,1365.036598,1352.333909,1355.242281,1347.182469,1349.681355,1347.173564,1338.731432,1342.823136,1342.767603,1337.76717,1331.681604,1333.078451,1334.301076,1331.570682,1325.092769,1327.707443,1311.448169,1318.968972,1317.336886,1313.694239,1297.235888,1297.013003,1289.323797,1289.885826,1283.736624,1277.455875,1277.521085,1269.173952,1267.788059,1263.577538,1260.326998,1254.803074,1269.809159,1266.117363, };
    double atemp[2 * (rank_ + 1)] = {0}, b[rank_ + 1] = {0}, a[rank_ + 1][rank_ + 1];
    int i, j, k;

    printf("%ld %ld\n", sizeof(x) / sizeof(double), sizeof(y) / sizeof(double));
    for(i = 0; i < maxn; i++){  //
        atemp[1] += x[i];
        atemp[2] += x[i]*x[i];
        atemp[3] += x[i]*x[i]*x[i];
        atemp[4] += x[i]*x[i]*x[i]*x[i];
        atemp[5] += x[i]*x[i]*x[i]*x[i]*x[i];
        atemp[6] += x[i]*x[i]*x[i]*x[i]*x[i]*x[i];
        b[0] += y[i];
        b[1] += x[i] * y[i];
        b[2] += x[i]*x[i]* y[i];
        b[3] += x[i]*x[i]*x[i] * y[i];
    }

    atemp[0] = maxn;
    /*
    for(i = 0; i <= 2 * rank_; i++)  printf("atemp[%d] = %f\n", i, atemp[i]);
    printf("\n");
    for(i = 0; i <= rank_; i++)  printf("b[%d] = %f\n", i, b[i]);
    printf("\n");
    */
    for(i = 0; i < rank_ + 1; i++){  //构建线性方程组系数矩阵，b[]不变
        k = i;
        for(j = 0; j < rank_ + 1; j++)  a[i][j] = atemp[k++];
    }
    /*
    for(i = 0; i < rank_ + 1; i++){
        for(j = 0; j < rank_ + 1; j++)  printf("a[%d][%d] = %-17f  ", i, j, a[i][j]);
        printf("\n");
    }
    printf("\n");
    */

    //以下为高斯列主元消去法解线性方程组
    for(k = 0; k < rank_ + 1 - 1; k++){  //n - 1列
        int column = k;
        double mainelement = a[k][k];

        for(i = k; i < rank_ + 1; i++)  //找主元素
            if(fabs(a[i][k]) > mainelement){
                mainelement = fabs(a[i][k]);
                column = i;
            }
        for(j = k; j < rank_ + 1; j++){  //交换两行
            double atemp = a[k][j];
            a[k][j] = a[column][j];
            a[column][j] = atemp;
        }
        double btemp = b[k];
        b[k] = b[column];
        b[column] = btemp;

        for(i = k + 1; i < rank_ + 1; i++){  //消元过程
            double Mik = a[i][k] / a[k][k];
            for(j = k; j < rank_ + 1; j++)  a[i][j] -= Mik * a[k][j];
            b[i] -= Mik * b[k];
        }
    }
    /*
    for(i = 0; i < rank_ + 1; i++){  //经列主元高斯消去法得到的上三角阵(最后一列为常系数)
        for(j = 0; j < rank_ + 1; j++)  printf("%20f", a[i][j]);
        printf("%20f\n", b[i]);
    }
    printf("\n");
    */
    b[rank_ + 1 - 1] /= a[rank_ + 1 - 1][rank_ + 1 - 1];  //回代过程
    for(i = rank_ + 1 - 2; i >= 0; i--){
        double sum = 0;
        for(j = i + 1; j < rank_ + 1; j++)  sum += a[i][j] * b[j];
        b[i] = (b[i] - sum) / a[i][i];
    }
    //高斯列主元消去法结束

    printf("P(x) = %f%+fx%+fx^2%+fx^3\n\n", b[0], b[1], b[2], b[3]);
    /*
    for(i = 0; i < maxn; i++){  //误差比较
        double temp = b[0] + b[1] * x[i] + b[2] * x[i] * x[i] + b[3] * x[i] * x[i] * x[i];
        printf("%f    %f    error: %f\n", y[i], temp, temp - y[i]);
    }
    */
    double res = Calc(-600, 600, b[0], b[1], b[2], b[3]);
    printf("res=%lf\n", res);
    return 0;
}